bulk-org
Raphael Kubler
12/15/2023
Introduction
Integrated bulk RNA-seq analysis (MSN9, MSN38, WTC11) of cerebral organoid samples in conditions co-cultured (with/without microglia) and stimulation (ctr/LPS/IFNy).
V1.83: corrected metadata and sample gene expression so that the analysis uses the right samples with proper annotation. V1.84: corrected the gene set used for PC variance analysis across groups.
CO vs COiMG
Organoids co-cultured with vs without microglia.
Variance analysis
PCAplot(pca_cor3_mod1, "COiMg", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='COiMG',y='Phagocytosis', fill='COiMG')) +
stat_compare_means(method = "t.test", label.y = 10, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Phagocytosis differential PC expression",x="COiMG", y = "PC1 Phagocytosis") +
theme_classic()
Microglia gene heatmap in ctrl samples
#for control only
ctr.meta <- meta[meta$condition %in% 'ctrl',]
ctr.meta$COiMg[ctr.meta$COiMg %in% 'yes'] <- 'COiMg'
ctr.meta$COiMg[ctr.meta$COiMg %in% 'no'] <- 'CO'
#horizontal
ra <- rowAnnotation(
COiMG = ctr.meta$COiMg[ctr.meta$condition %in% 'ctrl'],
col = list(
COiMG = c("COiMg"='tomato','CO'='lightblue')),
show_annotation_name = T,
show_legend = T)
ctr <- t(batch.rem3_mod1[,meta$condition %in% 'ctrl'])
hm_counts <- scale(ctr[,colnames(ctr) %in% gene_panels_subset$Microglia])
ComplexHeatmap::Heatmap(hm_counts,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(hm_counts), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
DEA
Volcano
volcano_plot <- function(res, title = NULL, subtitle = NULL, annotate_by = NULL, type ='ALS', ymax1 = 7.5, ymax2 = 8, xmax1 = -4.5, xmax2 = 4.5){
res <-
mutate(res,
sig = case_when(
padj >= 0.05 ~ "non_sig",
padj < 0.05 & abs(log2FoldChange) < 1 ~ "sig",
padj < 0.05 & abs(log2FoldChange) >= 1 ~ "sig - strong"
)) %>%
mutate(direction = ifelse(log2FoldChange > 0, "up", "down")) %>%
mutate(log2FoldChange = case_when(
log2FoldChange > 10 ~ Inf,
log2FoldChange < -10 ~ -Inf,
TRUE ~ log2FoldChange
)) %>%
mutate(class = paste(sig, direction))
if( type == "ALS"){
xpos <- 0.5
ymax <- ymax1
xlim <- c(xmax1,xmax2)
}else{
xpos <- 0.025
ymax <- ymax2
xlim <- c(-0.042, 0.042)
}
de_tally <- group_by(res, sig, direction, class) %>% tally() %>%
filter(sig != "non_sig") %>%
mutate( position = ifelse(sig == "sig", xpos, 2) ) %>%
mutate(position = ifelse( direction == "down", -1 * position, position)) %>%
mutate(n = formatC(n, format="f", big.mark=",", digits=0))
plot <- res %>%
mutate( pvalue = ifelse( pvalue < 1e-90, Inf, pvalue)) %>% #threshold at 1e16
ggplot(aes(x = log2FoldChange, y = -log10(pvalue))) +
#geom_point(aes(colour = class ), size = 0.5) +
rasterise(geom_point(aes(colour = class ), size = 0.5), dpi = 300) +
scale_colour_manual(values = c("non_sig up" = "gray",
"non_sig down" = "gray",
"sig up" = "#EB7F56",
"sig - strong up" = "#B61927",
"sig down" = "#4F8FC4",
"sig - strong down" = "dodgerblue4"
)) +
theme_bw() +
labs(y = expression(-log[10]~P~value), x = expression(log[2]~"(fold change)"), title = title, subtitle = subtitle) +
guides(colour = FALSE) +
scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.border = element_blank(),
axis.ticks = element_line(colour = "black")
) +
geom_text(fontface = "bold", data = de_tally, aes(x = position, y = ymax - 0.5, label = n, colour = class), size = 4) +
scale_x_continuous(limits = xlim)
if(!is.null(annotate_by)){
plot <- plot +
ggrepel::geom_text_repel(
fontface = "italic",
data = filter(res, symbol %in% annotate_by),
aes(x = log2FoldChange, y = -log10(pvalue), label = symbol),
min.segment.length = unit(0, "lines"),
size = 2.3) +
geom_point(
data = filter(res, symbol %in% annotate_by), size = 0.8, colour = "black"
) +
geom_point(aes(colour = class ),
data = filter(res, symbol %in% annotate_by), size = 0.6
)
}
return(plot)
}
results$COiMG$symbol <- rownames(results$COiMG)
volcano_plot(data.frame(results$COiMG), title = 'CO vs COiMg',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Top 50 DEGs heatmap
top50.coimg <- rbind(results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names <- rownames(top50.coimg)
ra2 <- rowAnnotation(
COiMG = meta$COiMg,
col = list(
COiMG = c("yes"='tomato','no'='lightblue')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra2,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
coimg.deg <- rbind(results$COiMG[rownames(results$COiMG) %in% na.omit(all_res$COiMG$DEGs[,1]),],
results$COiMG[rownames(results$COiMG) %in% na.omit(all_res$COiMG$DEGs[,2]),])
df <- cbind('DEG'=rownames(coimg.deg),'lFC'=round(coimg.deg$log2FoldChange,2),'Cilia-associated'=rownames(coimg.deg) %in% organoid$CBC)
createDT(data.frame(df), "Status", scrollY=1000)GO
dotplot(all_res$COiMG$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$COiMG$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
cgenes <- readxl::read_xlsx('Gene lists for organoid paper.xlsx', col_names = T)
m <- data.frame(t(toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`))[toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`)) %in% rownames(bulkorg_counts2)]), check.names = F)
for(i in names(cgenes)[-1]){
k <- na.omit(cgenes[i])
s <- k[as.matrix(k)[,1] %in% rownames(bulkorg_counts2),]
m <- rbind.fill(data.frame(m),data.frame(t(s)))
}
rownames(m) <- names(cgenes)
m <- t(m)
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
IFNy stimulation
Let’s do this for IFNy stimulation
Variance analysis
PCAplot(pca_cor3_mod1, "condition", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
DEA
Volcano
results$IFNg$symbol <- rownames(results$IFNg)
volcano_plot(data.frame(results$IFNg), title = 'CTR vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Top 50 DEGs heatmap
top50.ifng <- rbind(results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names2 <- rownames(top50.ifng)
ra3 <- rowAnnotation(
Stimulation = meta$condition,
col = list(
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra3,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
IFNg.deg <- rbind(results$IFNg[rownames(results$IFNg) %in% na.omit(all_res$IFNg$DEGs[,1]),],
results$IFNg[rownames(results$IFNg) %in% na.omit(all_res$IFNg$DEGs[,2]),])
df <- cbind('DEG'=rownames(IFNg.deg),'lFC'=round(IFNg.deg$log2FoldChange,2),'Cilia-associated'=rownames(IFNg.deg) %in% organoid$CBC)
createDT(data.frame(df), "Status", scrollY=1000)GO
dotplot(all_res$IFNg$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$IFNg$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
For organoid cell-types.
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For NDD gene lists.
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
DEGs$`Down-regulated`[DEGs$`Down-regulated` %in% as.character(na.omit(m[,12]))]## [1] "CYP7B1"For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
Cell-type marker analysis
Cell-types with only 1 gene found back in our expression data carry principal component expression = 0.
Heatmaps of gene & PC expression
#sign reversal
ctr.meta <- meta[meta$condition %in% 'ctrl',]
ctr.expr <- batch.rem3_mod1[,colnames(batch.rem3_mod1) %in% rownames(ctr.meta)]
organoid2 <- t(rbind.fill(data.frame(t(organoid)),data.frame(t(gene_panels_subset$Microglia))))
colnames(organoid2) <- c(colnames(organoid),'Microglia')
PC.covs <- list()
for(i in colnames(organoid2)){
pca.cov <- prcomp(t(ctr.expr[rownames(ctr.expr) %in% na.omit(organoid2[,colnames(organoid2) %in% i]),]))
PC.covs[[i]] <- -pca.cov$x
#if (table(rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1])[[2]] == 1){
# PC.covs[,i] <- t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],])
#}
}
for (i in names(PC.covs)){
if (length(colnames(PC.covs[[i]])) > 1){
df1 <- reshape::melt(PC.covs[[i]])
colnames(df1) <- c('Sample','PC','Rotation')
df1$COiMg <- ctr.meta$COiMg
p1 <- ggplot(df1,aes_string(x='PC',y='Rotation', fill='COiMg')) +
#stat_compare_means(method = "t.test", label.y = 6) +
geom_boxplot() +
labs(title=paste0('PC expression of ', i, ' markers'),x="Principal component", y = "PC expression") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
plot(p1)
}
}
Boxplots
Let’s make boxplots instead
#add variance analysis of cell-type markersDissecting DEA results
CO-IFNy
These are the results of the DEA of CTR vs IFNy in CO (not co-cultured with microglia) only. Cilia-associated genes remain downregulated. We find this as well in CO vs COiMG and more so also in our single-cell data, which did not have any IFNy-stimulated samples. Meaning: this is a result robust to IFNy AND microglia co-culture separately.
df <- cbind('DEG'=rownames(IFNy.res),'lFC'=round(IFNy.res$log2FoldChange,2),'Cilia-associated'=rownames(IFNy.res) %in% organoid$CBC)
createDT(data.frame(df), "Stimulation", scrollY=1000)Volcano
new_IFNy$symbol <- rownames(new_IFNy)
volcano_plot(data.frame(new_IFNy), title = 'ctrl vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
CO-COiMg
This is to show that the DEA effect of cilia-associated genes remains present in our bulk-seq analysis when we subset to control samples only and test CO vs COiMg.
df <- cbind('DEG'=rownames(coimg.res),'lFC'=round(coimg.res$log2FoldChange,2),'Cilia-associated'=rownames(coimg.res) %in% organoid$CBC)
createDT(data.frame(df), "COiMG", scrollY=1000)Volcano
new_coimg$symbol <- rownames(new_coimg)
volcano_plot(data.frame(new_coimg), title = 'CO vs COiMg',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Correlation of lFCs
overlap <- rownames(coimg.res[rownames(coimg.res) %in% rownames(IFNy.res),])
cor(coimg.res[rownames(coimg.res) %in% overlap,]$log2FoldChange,IFNy.res[rownames(IFNy.res) %in% overlap,]$log2FoldChange)## [1] 0.8821443createDT(cbind('Genes'=overlap,
'COiMg'=coimg.res[rownames(coimg.res) %in% overlap,]$log2FoldChange,
'IFNy'=IFNy.res[rownames(IFNy.res) %in% overlap,]$log2FoldChange), "COiMG", scrollY=1000)By Line analysis
Correlation of DEA results across lines by condition
IFNy
'IFNy_overlap_wtc11Xmsn9:'## [1] "IFNy_overlap_wtc11Xmsn9:"line.res$DEG_IFNy$overlap_wtc11Xmsn9## [1] "BCAN" "HPGD" "CFH" "LAP3" "CASP10" "CD38"
## [7] "NOS2" "MMP25" "IL32" "ETV7" "MRC2" "ALOX5"
## [13] "MVP" "NUB1" "CD74" "TYMP" "FAS" "CD44"
## [19] "BTN3A1" "VRK2" "TNFRSF1B" "PI4K2B" "TNC" "LCP2"
## [25] "GUCA1A" "EIF2AK2" "PRDM1" "PARP12" "OGFR" "CASP8"
## [31] "SBNO2" "MTHFD2" "SP100" "IFI35" "BCL3" "RAB27A"
## [37] "TESK2" "TBX21" "ARHGAP15" "IL4R" "LAMP3" "EDN1"
## [43] "SP140" "CEACAM1" "FYB1" "WDFY1" "CETP" "EBF4"
## [49] "OAS1" "ICAM1" "FLT3LG" "PSME1" "FKBP5" "JAK2"
## [55] "IL12RB1" "HSD3B7" "APOL4" "APOL1" "LGMN" "PCK2"
## [61] "PSME2" "CD40" "TRIB3" "TLDC2" "SAMHD1" "SMCHD1"
## [67] "PRPS2" "TLR8" "ELF4" "HTR2A" "TNFSF13B" "FOXF1"
## [73] "RIPK2" "NBN" "GSDMD" "SLC1A5" "SIGLEC8" "PLEKHA4"
## [79] "NAMPT" "ZC3HAV1" "MET" "NOD1" "GIMAP2" "TRIM14"
## [85] "RIGI" "ACTA2" "MAP3K8" "SORCS1" "LGALS3BP" "DHX58"
## [91] "FAM20A" "CPZ" "UNC93B1" "SLC15A3" "TCIRG1" "CLEC2B"
## [97] "OAS3" "OAS2" "ADGRD1" "BTN3A3" "TRIM38" "TENT5A"
## [103] "ITK" "SLC12A7" "BCL6" "XRN1" "GNB4" "CISH"
## [109] "OTOF" "IFIH1" "STAT1" "FBXO6" "GBP3" "GBP1"
## [115] "IFIT3" "IFIT2" "CD274" "TRIM25" "ABCC11" "BCL2L14"
## [121] "TNFSF10" "GLIPR2" "SSPN" "ADGRE5" "NMI" "BATF3"
## [127] "DNPEP" "ZNFX1" "RNF114" "ZBP1" "MT2A" "IRF1"
## [133] "RBCK1" "IGFLR1" "IFI6" "L3HYPDH" "HELB" "OR7C1"
## [139] "STEAP4" "ZFP36" "ADM2" "APOL3" "APOL2" "APOBEC3F"
## [145] "ACKR4" "AIPL1" "RHBDF2" "STARD8" "PLVAP" "BST2"
## [151] "SAG" "HELZ2" "SHFL" "CCL25" "IDO1" "RFTN1"
## [157] "CHSY1" "GCH1" "PODNL1" "TRIM21" "TRIM22" "XAF1"
## [163] "EPSTI1" "PLAAT4" "C1QTNF6" "PSRC1" "NGF" "RSAD2"
## [169] "IL15RA" "ANXA1" "OASL" "PRPH" "GLUL" "SP110"
## [175] "DOCK10" "PHF11" "RTP4" "KLF4" "IL18BP" "DDX60"
## [181] "CASP1" "SQOR" "IFI44L" "IFI44" "STAMBPL1" "PARP9"
## [187] "HERC6" "CXCL9" "TAPBPL" "FBLN5" "WARS1" "PML"
## [193] "HAPLN3" "MARVELD3" "NLRC5" "IRF8" "ZMYND15" "SECTM1"
## [199] "PMAIP1" "IFITM3" "MOB3C" "FCGR2A" "ADAMTSL4" "RALB"
## [205] "C1QL2" "PLAC8" "CYP4V2" "OSMR" "TNFRSF21" "TMEM140"
## [211] "SERPING1" "KCTD14" "FLI1" "AP1S3" "PSTPIP2" "ANKRD22"
## [217] "IFIT5" "GBP5" "L3MBTL4" "TTC39B" "GPR78" "MOV10"
## [223] "PIK3AP1" "CLIC2" "ART3" "MAP3K7CL" "UBE2L6" "MX1"
## [229] "TNFRSF14" "MEGF11" "SLAMF8" "IFNAR2" "CBR3" "C1R"
## [235] "PTH1R" "LY6E" "MFSD12" "CXCL16" "SCIMP" "PDPN"
## [241] "GBP2" "GBP4" "ATF3" "PIGR" "CTSS" "S100A11"
## [247] "FZD5" "NUAK2" "IFI16" "CDCP1" "CCR1" "DTX3L"
## [253] "SHISA5" "MST1R" "TMEM144" "IL15" "ELOVL7" "F2RL2"
## [259] "ERAP1" "ERAP2" "TLR3" "DACT2" "CASP7" "IFI27"
## [265] "C2" "PLEKHF1" "TVP23A" "B2M" "SLFN5" "FAM111A"
## [271] "TBC1D2B" "NOD2" "PRRT2" "BATF2" "IRF2" "XIRP1"
## [277] "FILIP1L" "TAP1" "MLKL" "RHOH" "LGI3" "STAT3"
## [283] "SNTB2" "VAMP5" "LGALS9" "PXDC1" "CXCL10" "CXCL11"
## [289] "PLEKHA2" "PROKR1" "TRIM56" "TCAF2" "STAT2" "GPR82"
## [295] "TMEM51" "ADRA1D" "ISG20" "MBOAT1" "CEBPB" "SLFN11"
## [301] "MYD88" "PARP14" "SLC2A14" "CD7" "RNF213" "MARCHF3"
## [307] "DES" "GRAMD2A" "A2M" "NUPR1" "FBXO39" "SAMD9L"
## [313] "ODF3B" "PARP10" "EXOC3L1" "CIITA" "GPBAR1" "PCGF5"
## [319] "DDX60L" "UBA7" "C1S" "LHFPL1" "GBP6" "MX2"
## [325] "LHFPL6" "ASCL2" "CSF1" "SOCS3" "STING1" "PRSS55"
## [331] "RBM43" "USP18" "MAFF" "TNFAIP2" "SOCS1" "SP140L"
## [337] "MUC1" "IRF7" "IFIT1" "TRIM69" "IFITM1" "ZNF267"
## [343] "GLDN" "BTN3A2" "ISG15" "RUFY4" "PLSCR1" "DPYD"
## [349] "CALHM6" "INSYN2A" "PCDH18" "CXCL17" "HLA-DRB1" "CD47"
## [355] "OR52K1" "ARID5A" "CASP4" "ACSL5" "GSTK1" "MAP1LC3C"
## [361] "TEAD4" "TGM2" "HLA-DOA" "HLA-DMA" "PSMB8" "TAP2"
## [367] "HLA-DRA" "CARD16" "HLA-C" "HLA-E" "HLA-F" "PSMB10"
## [373] "SAMD9" "ATP10A" "COL6A6" "HLA-A" "GIMAP1" "SLC35F6"
## [379] "UBD" "IRF9" "SMTNL1" "IFI30" "CEBPD" "APOL6"
## [385] "HLA-DPB1" "C4B" "PATL2" "HLA-DPA1" "CYP21A2" "TAPBP"
## [391] "HLA-B" "OR2I1P" "APOBEC3G" "PSMB9" "HLA-DOB" "HLA-DMB"
## [397] "CFB" "APOBEC3D" "APOBEC3C" "C4A" "PRSS51" "TRIL"
## [403] "CTSO" "ZNF878" "TRIM34" "CCDC194" "CCL5" "PAGR1"
## [409] "TCAF2C"print('Correlation of IFNy results between WTC11 & MSN9')## [1] "Correlation of IFNy results between WTC11 & MSN9"cor(line.res$DEG_IFNy[['WTC11']][rownames(line.res$DEG_IFNy[['WTC11']]) %in% line.res$DEG_IFNy$overlap_wtc11Xmsn9,]$log2FoldChange,
line.res$DEG_IFNy[['MSN9']][rownames(line.res$DEG_IFNy[['MSN9']]) %in% line.res$DEG_IFNy$overlap_wtc11Xmsn9,]$log2FoldChange)## [1] 0.8980987'###################'## [1] "###################"'IFNy_overlap_msn9Xmsn38:'## [1] "IFNy_overlap_msn9Xmsn38:"line.res$DEG_IFNy$overlap_msn9Xmsn38## [1] "CFTR" "PCDH11X" "HPGD" "CFH" "LAP3" "CASP10"
## [7] "CD38" "NOS2" "MMP25" "IL32" "ETV7" "MRC2"
## [13] "MVP" "NUB1" "CD74" "HGF" "TYMP" "FAS"
## [19] "CD44" "BTN3A1" "TNFRSF1B" "TNC" "JADE2" "LCP2"
## [25] "GUCA1A" "EIF2AK2" "PRDM1" "PARP12" "OGFR" "CASP8"
## [31] "SBNO2" "MTHFD2" "SP100" "IFI35" "BCL3" "ALPK1"
## [37] "TBX21" "ARHGAP15" "LAMP3" "SP140" "CEACAM1" "FYB1"
## [43] "WDFY1" "CETP" "OAS1" "ICAM1" "FLT3LG" "PSME1"
## [49] "JAK2" "HSD3B7" "APOL4" "APOL1" "PCK2" "PSME2"
## [55] "CD40" "TRIB3" "TLDC2" "SAMHD1" "SMCHD1" "TNFSF13B"
## [61] "LAPTM4B" "IL7" "GSDMD" "RELB" "ZC3HAV1" "GIMAP2"
## [67] "TRIM14" "RIGI" "ACTA2" "LGALS3BP" "DHX58" "FAM20A"
## [73] "CPZ" "SLC15A3" "CLEC2B" "OAS3" "OAS2" "BTN3A3"
## [79] "TRIM38" "TENT5A" "ITK" "SLC12A7" "BCL6" "CISH"
## [85] "OTOF" "IFIH1" "STAT1" "FBXO6" "GBP3" "GBP1"
## [91] "IFIT3" "IFIT2" "CD274" "TRIM25" "ABCC11" "TNFSF10"
## [97] "NMI" "BATF3" "DNPEP" "ZNFX1" "MT2A" "IRF1"
## [103] "RBCK1" "IGFLR1" "STAT5A" "IFI6" "HELB" "STEAP4"
## [109] "ZFP36" "APOL3" "APOL2" "APOBEC3F" "STARD8" "PLVAP"
## [115] "BST2" "SAG" "HELZ2" "SHFL" "CCL25" "IDO1"
## [121] "GCH1" "TRIM21" "TRIM5" "TRIM22" "XAF1" "EPSTI1"
## [127] "PLAAT4" "GIMAP4" "PSRC1" "RSAD2" "CMPK2" "IL15RA"
## [133] "OASL" "PRPH" "STX11" "SP110" "DOCK10" "PHF11"
## [139] "RTP4" "KLF4" "IL18BP" "DDX60" "CASP1" "SQOR"
## [145] "IFI44L" "IFI44" "STAMBPL1" "PARP9" "HERC6" "HERC5"
## [151] "CXCL9" "ETV6" "TAPBPL" "RAB20" "FBLN5" "WARS1"
## [157] "PML" "HAPLN3" "MARVELD3" "NLRC5" "IRF8" "SECTM1"
## [163] "PMAIP1" "IFITM3" "MOB3C" "FCGR2A" "TIFA" "CYP4V2"
## [169] "OSMR" "SERPING1" "AP1S3" "ANKRD22" "IFIT5" "GBP5"
## [175] "L3MBTL4" "MOV10" "PIK3AP1" "CLIC2" "ART3" "MAP3K7CL"
## [181] "UBE2L6" "KCNJ15" "MX1" "TNFRSF14" "MEGF11" "MITD1"
## [187] "SLAMF8" "IFNAR2" "RUNX1" "CBR3" "C1R" "ICOSLG"
## [193] "LY6E" "CXCL16" "PDPN" "GBP2" "GBP4" "ATF3"
## [199] "NEURL3" "CTSS" "ARHGAP25" "FZD5" "IFI16" "CDCP1"
## [205] "DTX3L" "TMEM144" "IL15" "ERAP1" "ERAP2" "TLR3"
## [211] "KCNV1" "CASP7" "IFI27" "C2" "PLEKHF1" "TVP23A"
## [217] "B2M" "TBC1D2B" "NOD2" "PRRT2" "BATF2" "IRF2"
## [223] "FILIP1L" "TAP1" "MLKL" "STAT3" "VAMP5" "LGALS9"
## [229] "SIGLEC7" "RGS14" "CXCL10" "CXCL11" "COL22A1" "PLEKHA2"
## [235] "PROKR1" "TRIM56" "TCAF2" "STAT2" "TMEM51" "ISG20"
## [241] "CEBPB" "SLFN11" "MYD88" "PARP14" "CD7" "RNF213"
## [247] "MARCHF3" "DES" "RMI2" "A2M" "FBXO39" "SAMD9L"
## [253] "PARP10" "EXOC3L1" "PXT1" "CIITA" "DDX60L" "UBA7"
## [259] "C1S" "GBP6" "MX2" "ASCL2" "CSF1" "SOCS3"
## [265] "USP18" "MAFF" "TNFAIP2" "SOCS1" "SP140L" "IRF7"
## [271] "IFIT1" "TRIM69" "IFITM1" "BTN3A2" "ISG15" "RUFY4"
## [277] "PLSCR1" "CALHM6" "INSYN2A" "HLA-DRB1" "ARID5A" "CASP4"
## [283] "ACSL5" "GSTK1" "MAP1LC3C" "TEAD4" "HLA-DRB5" "TGM2"
## [289] "HLA-DOA" "HLA-DMA" "PSMB8" "TAP2" "HLA-DRA" "CARD16"
## [295] "MICB" "HLA-C" "HLA-E" "HLA-F" "PSMB10" "SAMD9"
## [301] "ATP10A" "HLA-A" "UBD" "IRF9" "SMTNL1" "IFI30"
## [307] "CEBPD" "APOL6" "HLA-DPB1" "C4B" "PATL2" "HLA-DPA1"
## [313] "CYP21A2" "TAPBP" "HLA-B" "ANKRD65" "OR2I1P" "APOBEC3G"
## [319] "PSMB9" "HLA-DOB" "HLA-DMB" "CFB" "APOBEC3D" "APOBEC3C"
## [325] "C4A" "PRSS51" "TRIL" "ZNF878" "TRIM34" "CCDC194"
## [331] "CCL5"print('Correlation of IFNy results between MSN38 & MSN9')## [1] "Correlation of IFNy results between MSN38 & MSN9"cor(line.res$DEG_IFNy[['MSN38']][rownames(line.res$DEG_IFNy[['MSN38']]) %in% line.res$DEG_IFNy$overlap_msn9Xmsn38,]$log2FoldChange,
line.res$DEG_IFNy[['MSN9']][rownames(line.res$DEG_IFNy[['MSN9']]) %in% line.res$DEG_IFNy$overlap_msn9Xmsn38,]$log2FoldChange)## [1] 0.8857025'###################'## [1] "###################"'IFNy_overlap_wtc11Xmsn38:'## [1] "IFNy_overlap_wtc11Xmsn38:"line.res$DEG_IFNy$overlap_msn38Xwtc11## [1] "HPGD" "GRIN2A" "ERVMER34-1" "CFH" "LAP3"
## [6] "CASP10" "CD38" "NOS2" "MMP25" "IL32"
## [11] "ETV7" "MRC2" "MVP" "NUB1" "WAS"
## [16] "CD74" "SAMD4A" "TYMP" "FAS" "CD44"
## [21] "BTN3A1" "TNFRSF1B" "TNC" "LCP2" "ARAP2"
## [26] "GUCA1A" "EIF2AK2" "PRDM1" "PARP12" "OGFR"
## [31] "LZTS1" "CASP8" "SBNO2" "MTHFD2" "SP100"
## [36] "IFI35" "BCL3" "ASNS" "TBX21" "CACNG5"
## [41] "ARHGAP15" "LAMP3" "SP140" "CEACAM1" "IL12RB2"
## [46] "FYB1" "WDFY1" "CETP" "OAS1" "ICAM1"
## [51] "FLT3LG" "PSME1" "RNF31" "JAK2" "HSD3B7"
## [56] "APOL4" "APOL1" "PCK2" "PSME2" "CD40"
## [61] "TRIB3" "TLDC2" "SAMHD1" "SMCHD1" "TNFSF13B"
## [66] "SLC7A5" "GSDMD" "ZC3HAV1" "CAV2" "CAV1"
## [71] "GIMAP2" "TRIM14" "RIGI" "ACTA2" "LGALS3BP"
## [76] "DHX58" "FAM20A" "CPZ" "SLC15A3" "CARS1"
## [81] "CLEC2B" "OAS3" "OAS2" "BTN3A3" "NCOA7"
## [86] "TRIM38" "TENT5A" "ITK" "SLC12A7" "BCL6"
## [91] "CISH" "OTOF" "IFIH1" "STAT1" "RNF19B"
## [96] "FBXO6" "GBP3" "GBP1" "IFIT3" "IFIT2"
## [101] "CD274" "PTK2B" "TRIM25" "ABCC11" "TNFSF10"
## [106] "NRK" "NMI" "BATF3" "DNPEP" "ZNFX1"
## [111] "MT2A" "IRF1" "C3" "RBCK1" "IGFLR1"
## [116] "IFI6" "HELB" "FGL2" "STEAP4" "ZFP36"
## [121] "APOL3" "APOL2" "APOBEC3F" "CHAC1" "STARD8"
## [126] "PLVAP" "BST2" "SAG" "HELZ2" "SHFL"
## [131] "CCL25" "IDO1" "GCH1" "TRIM21" "TRIM22"
## [136] "XAF1" "EPSTI1" "PLAAT4" "DGLUCY" "PSRC1"
## [141] "RSAD2" "IL15RA" "OASL" "PRRG4" "PRPH"
## [146] "SP110" "DOCK10" "ALDH1L2" "PHF11" "RTP4"
## [151] "KLF4" "RNF144B" "IL18BP" "DDX60" "CASP1"
## [156] "SQOR" "IFI44L" "IFI44" "STAMBPL1" "PARP9"
## [161] "HERC6" "CXCL9" "TAPBPL" "JDP2" "FBLN5"
## [166] "WARS1" "PML" "HAPLN3" "MARVELD3" "NLRC5"
## [171] "IRF8" "SECTM1" "PMAIP1" "IFITM3" "MOB3C"
## [176] "FCGR2A" "CYP4V2" "OSMR" "TACC1" "NACC2"
## [181] "SERPING1" "AP1S3" "RASGRP3" "ANKRD22" "IFIT5"
## [186] "LGI4" "GBP5" "L3MBTL4" "MOV10" "PIK3AP1"
## [191] "CLIC2" "ART3" "MAP3K7CL" "UBE2L6" "MX1"
## [196] "TNFRSF14" "MEGF11" "SLAMF8" "IFNAR2" "CBR3"
## [201] "C1R" "ADAR" "LY6E" "CXCL16" "PDPN"
## [206] "GBP2" "GBP4" "VCAM1" "ATF3" "CTSS"
## [211] "FZD5" "IFI16" "CDCP1" "DTX3L" "LIPH"
## [216] "TMEM144" "IL15" "ERAP1" "ERAP2" "TLR3"
## [221] "CASP7" "OTOGL" "IFI27" "C2" "PLEKHF1"
## [226] "TVP23A" "B2M" "TBC1D2B" "NOD2" "PRRT2"
## [231] "BATF2" "IRF2" "FILIP1L" "TAP1" "MLKL"
## [236] "STAT3" "VAMP5" "LGALS9" "ATF5" "CXCL10"
## [241] "CXCL11" "PLEKHA2" "PROKR1" "TRIM56" "TCAF2"
## [246] "STAT2" "TMEM51" "ISG20" "CEBPB" "SLFN11"
## [251] "MYD88" "PARP14" "PARP15" "C1QB" "C1QA"
## [256] "CD7" "RNF213" "MARCHF3" "DES" "DDIT3"
## [261] "A2M" "FBXO39" "SAMD9L" "PARP10" "EXOC3L1"
## [266] "CIITA" "DDX60L" "UBA7" "C1S" "GBP6"
## [271] "MX2" "ASCL2" "FAM162B" "CSF1" "SOCS3"
## [276] "USP18" "MAFF" "TNFAIP2" "SOCS1" "SP140L"
## [281] "IRF7" "IFIT1" "TRIM69" "IFITM1" "BTN3A2"
## [286] "ISG15" "EIF4EBP1" "RUFY4" "PLSCR1" "CALHM6"
## [291] "INSYN2A" "ANKRD34B" "HLA-DRB1" "ARID5A" "CASP4"
## [296] "ACSL5" "GSTK1" "MAP1LC3C" "TEAD4" "ALPK2"
## [301] "TGM2" "HLA-DOA" "HLA-DMA" "PSMB8" "TAP2"
## [306] "HLA-DRA" "CARD16" "LST1" "HLA-C" "HLA-E"
## [311] "HLA-F" "PSMB10" "SAMD9" "ATP10A" "HLA-A"
## [316] "UBD" "IRF9" "SMTNL1" "IFI30" "CEBPD"
## [321] "APOL6" "HLA-DPB1" "C4B" "PATL2" "HLA-DPA1"
## [326] "CYP21A2" "TAPBP" "HLA-B" "KIAA0040" "OR2I1P"
## [331] "APOBEC3G" "PSMB9" "HLA-DOB" "HLA-DMB" "CFB"
## [336] "APOBEC3D" "APOBEC3C" "C4A" "TMEM158" "PRSS51"
## [341] "TRIL" "ZNF878" "TRIM34" "CCDC194" "CCL5"
## [346] "PRAG1"print('Correlation of IFNy results between MSN38 & WTC11')## [1] "Correlation of IFNy results between MSN38 & WTC11"cor(line.res$DEG_IFNy[['MSN38']][rownames(line.res$DEG_IFNy[['MSN38']]) %in% line.res$DEG_IFNy$overlap_msn38Xwtc11,]$log2FoldChange,
line.res$DEG_IFNy[['WTC11']][rownames(line.res$DEG_IFNy[['WTC11']]) %in% line.res$DEG_IFNy$overlap_msn38Xwtc11,]$log2FoldChange)## [1] 0.904527LinexCondition
for (i in c('FOXJ1',
'CLDN5',
'IDO1',
'GBP5',
'MMP10',
'SOX2',
'RBFOX3')){
p <- ggplot(pca.df_upgraded,aes_string(x='LinexCondition',y=i, fill='LinexCondition')) +
#stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title=paste0("LinexCondition ", i, " expression"),x="LinexCondition", y = paste0(i,' expression')) +
theme_classic()
plot(p)
}
COiMg
'IFNy_overlap_wtc11Xmsn9:'## [1] "IFNy_overlap_wtc11Xmsn9:"line.res$DEG_COiMg$overlap_wtc11Xmsn9## [1] "BCAS1" "HSD17B2" "EPB41L4B" "RGCC" "KDR" "PDGFRA"
## [7] "WNT10A" "KCNH5" "THRB" "EFHB" "ATOH8" "TACR3"
## [13] "SOSTDC1" "PRIMA1" "PLD5" "AIFM3" "MFRP" "FGR"
## [19] "VCAN" "COL9A2" "FYB1" "RGS1" "HMOX1" "SIGLEC8"
## [25] "ODAM" "MDFI" "CD86" "TFCP2L1" "CYTIP" "PROC"
## [31] "PLEK" "CD48" "CCRL2" "SRGN" "PMEPA1" "FIBCD1"
## [37] "KANK4" "ADAMDEC1" "PTPN22" "FST" "HAVCR2" "SDS"
## [43] "TLR4" "CYP19A1" "IFI44" "GOLGA6D" "LEFTY2" "IL2RG"
## [49] "CDKN2B" "SLC34A2" "GOLGA6A" "CCKAR" "FBLN2" "TAGAP"
## [55] "PLD4" "CD52" "FPR1" "GPR34" "KLHL6" "LPL"
## [61] "MSC" "F2R" "UBA7" "AMTN" "HSH2D" "SERPINA1"
## [67] "DPP4" "SUCNR1" "TGM2" "AIF1" "ATP10A" "STK38L"
## [73] "GOLGA6B" "CLEC19A"print('Correlation of IFNy results between WTC11 & MSN9')## [1] "Correlation of IFNy results between WTC11 & MSN9"cor(line.res$DEG_COiMg[['WTC11']][rownames(line.res$DEG_COiMg[['WTC11']]) %in% line.res$DEG_COiMg$overlap_wtc11Xmsn9,]$log2FoldChange,
line.res$DEG_COiMg[['MSN9']][rownames(line.res$DEG_COiMg[['MSN9']]) %in% line.res$DEG_COiMg$overlap_wtc11Xmsn9,]$log2FoldChange)## [1] 0.5575205'###################'## [1] "###################"'IFNy_overlap_msn9Xmsn38:'## [1] "IFNy_overlap_msn9Xmsn38:"line.res$DEG_IFNy$overlap_msn9Xmsn38## [1] "CFTR" "PCDH11X" "HPGD" "CFH" "LAP3" "CASP10"
## [7] "CD38" "NOS2" "MMP25" "IL32" "ETV7" "MRC2"
## [13] "MVP" "NUB1" "CD74" "HGF" "TYMP" "FAS"
## [19] "CD44" "BTN3A1" "TNFRSF1B" "TNC" "JADE2" "LCP2"
## [25] "GUCA1A" "EIF2AK2" "PRDM1" "PARP12" "OGFR" "CASP8"
## [31] "SBNO2" "MTHFD2" "SP100" "IFI35" "BCL3" "ALPK1"
## [37] "TBX21" "ARHGAP15" "LAMP3" "SP140" "CEACAM1" "FYB1"
## [43] "WDFY1" "CETP" "OAS1" "ICAM1" "FLT3LG" "PSME1"
## [49] "JAK2" "HSD3B7" "APOL4" "APOL1" "PCK2" "PSME2"
## [55] "CD40" "TRIB3" "TLDC2" "SAMHD1" "SMCHD1" "TNFSF13B"
## [61] "LAPTM4B" "IL7" "GSDMD" "RELB" "ZC3HAV1" "GIMAP2"
## [67] "TRIM14" "RIGI" "ACTA2" "LGALS3BP" "DHX58" "FAM20A"
## [73] "CPZ" "SLC15A3" "CLEC2B" "OAS3" "OAS2" "BTN3A3"
## [79] "TRIM38" "TENT5A" "ITK" "SLC12A7" "BCL6" "CISH"
## [85] "OTOF" "IFIH1" "STAT1" "FBXO6" "GBP3" "GBP1"
## [91] "IFIT3" "IFIT2" "CD274" "TRIM25" "ABCC11" "TNFSF10"
## [97] "NMI" "BATF3" "DNPEP" "ZNFX1" "MT2A" "IRF1"
## [103] "RBCK1" "IGFLR1" "STAT5A" "IFI6" "HELB" "STEAP4"
## [109] "ZFP36" "APOL3" "APOL2" "APOBEC3F" "STARD8" "PLVAP"
## [115] "BST2" "SAG" "HELZ2" "SHFL" "CCL25" "IDO1"
## [121] "GCH1" "TRIM21" "TRIM5" "TRIM22" "XAF1" "EPSTI1"
## [127] "PLAAT4" "GIMAP4" "PSRC1" "RSAD2" "CMPK2" "IL15RA"
## [133] "OASL" "PRPH" "STX11" "SP110" "DOCK10" "PHF11"
## [139] "RTP4" "KLF4" "IL18BP" "DDX60" "CASP1" "SQOR"
## [145] "IFI44L" "IFI44" "STAMBPL1" "PARP9" "HERC6" "HERC5"
## [151] "CXCL9" "ETV6" "TAPBPL" "RAB20" "FBLN5" "WARS1"
## [157] "PML" "HAPLN3" "MARVELD3" "NLRC5" "IRF8" "SECTM1"
## [163] "PMAIP1" "IFITM3" "MOB3C" "FCGR2A" "TIFA" "CYP4V2"
## [169] "OSMR" "SERPING1" "AP1S3" "ANKRD22" "IFIT5" "GBP5"
## [175] "L3MBTL4" "MOV10" "PIK3AP1" "CLIC2" "ART3" "MAP3K7CL"
## [181] "UBE2L6" "KCNJ15" "MX1" "TNFRSF14" "MEGF11" "MITD1"
## [187] "SLAMF8" "IFNAR2" "RUNX1" "CBR3" "C1R" "ICOSLG"
## [193] "LY6E" "CXCL16" "PDPN" "GBP2" "GBP4" "ATF3"
## [199] "NEURL3" "CTSS" "ARHGAP25" "FZD5" "IFI16" "CDCP1"
## [205] "DTX3L" "TMEM144" "IL15" "ERAP1" "ERAP2" "TLR3"
## [211] "KCNV1" "CASP7" "IFI27" "C2" "PLEKHF1" "TVP23A"
## [217] "B2M" "TBC1D2B" "NOD2" "PRRT2" "BATF2" "IRF2"
## [223] "FILIP1L" "TAP1" "MLKL" "STAT3" "VAMP5" "LGALS9"
## [229] "SIGLEC7" "RGS14" "CXCL10" "CXCL11" "COL22A1" "PLEKHA2"
## [235] "PROKR1" "TRIM56" "TCAF2" "STAT2" "TMEM51" "ISG20"
## [241] "CEBPB" "SLFN11" "MYD88" "PARP14" "CD7" "RNF213"
## [247] "MARCHF3" "DES" "RMI2" "A2M" "FBXO39" "SAMD9L"
## [253] "PARP10" "EXOC3L1" "PXT1" "CIITA" "DDX60L" "UBA7"
## [259] "C1S" "GBP6" "MX2" "ASCL2" "CSF1" "SOCS3"
## [265] "USP18" "MAFF" "TNFAIP2" "SOCS1" "SP140L" "IRF7"
## [271] "IFIT1" "TRIM69" "IFITM1" "BTN3A2" "ISG15" "RUFY4"
## [277] "PLSCR1" "CALHM6" "INSYN2A" "HLA-DRB1" "ARID5A" "CASP4"
## [283] "ACSL5" "GSTK1" "MAP1LC3C" "TEAD4" "HLA-DRB5" "TGM2"
## [289] "HLA-DOA" "HLA-DMA" "PSMB8" "TAP2" "HLA-DRA" "CARD16"
## [295] "MICB" "HLA-C" "HLA-E" "HLA-F" "PSMB10" "SAMD9"
## [301] "ATP10A" "HLA-A" "UBD" "IRF9" "SMTNL1" "IFI30"
## [307] "CEBPD" "APOL6" "HLA-DPB1" "C4B" "PATL2" "HLA-DPA1"
## [313] "CYP21A2" "TAPBP" "HLA-B" "ANKRD65" "OR2I1P" "APOBEC3G"
## [319] "PSMB9" "HLA-DOB" "HLA-DMB" "CFB" "APOBEC3D" "APOBEC3C"
## [325] "C4A" "PRSS51" "TRIL" "ZNF878" "TRIM34" "CCDC194"
## [331] "CCL5"print('Correlation of IFNy results between MSN38 & MSN9')## [1] "Correlation of IFNy results between MSN38 & MSN9"cor(line.res$DEG_COiMg[['MSN38']][rownames(line.res$DEG_COiMg[['MSN38']]) %in% line.res$DEG_COiMg$overlap_msn9Xmsn38,]$log2FoldChange,
line.res$DEG_COiMg[['MSN9']][rownames(line.res$DEG_COiMg[['MSN9']]) %in% line.res$DEG_COiMg$overlap_msn9Xmsn38,]$log2FoldChange)## [1] 0.7178828'###################'## [1] "###################"'IFNy_overlap_msn9Xmsn38:'## [1] "IFNy_overlap_msn9Xmsn38:"line.res$DEG_COiMg$overlap_msn38Xwtc11## [1] "BCAS1" "DBX1" "FGR" "TMEM176A" "ITGAL" "TRAF3IP3"
## [7] "STAB1" "CD4" "BTK" "MRC2" "TYROBP" "ALOX5"
## [13] "SLC11A1" "RUNX3" "SLAMF7" "TNFRSF1B" "VCAN" "MSR1"
## [19] "CAPG" "ADAM28" "LCP2" "COL9A2" "ELN" "TBXAS1"
## [25] "GNA15" "CD84" "SPI1" "APBB1IP" "EDN1" "SP140"
## [31] "FYB1" "LAT2" "NOX4" "ADAMTS2" "SIGLEC1" "RGS1"
## [37] "LYZ" "UNC13D" "TREM2" "IL12RB1" "CYTH4" "HMOX1"
## [43] "NCF4" "CSF2RB" "HCK" "MXRA5" "ACP5" "PYCARD"
## [49] "IL4I1" "RASAL3" "SIGLEC8" "CD33" "PIK3CG" "TFEC"
## [55] "TMEM176B" "ENG" "ABI3" "COL1A1" "IL10RA" "SLC15A3"
## [61] "C11orf21" "BIN2" "ARHGDIB" "PTPN6" "MDFI" "LOX"
## [67] "CD86" "CYTIP" "ITGB6" "PROC" "PLEK" "NCF2"
## [73] "CD48" "PADI2" "SPP1" "CSF3R" "TGFBI" "CCRL2"
## [79] "TMIGD3" "SASH3" "SRGN" "BHLHE41" "ARHGAP9" "NCKAP1L"
## [85] "PMEPA1" "C3" "EVI2A" "ADGRE2" "FGL2" "RAC2"
## [91] "IRF5" "WDFY4" "CD68" "SIGLEC9" "APOC1" "LSP1"
## [97] "COL5A1" "FIBCD1" "CARD6" "ALOX5AP" "CHI3L1" "CHIT1"
## [103] "PRAM1" "GIMAP6" "GIMAP4" "ADAMDEC1" "CD180" "PTPN22"
## [109] "FST" "DOCK2" "SLC37A2" "HAVCR2" "SDS" "CD36"
## [115] "NT5E" "DYSF" "NPL" "LCP1" "IL1RN" "TLR4"
## [121] "SLCO2B1" "CASP1" "ITGA11" "PLCB2" "PARVG" "ACVRL1"
## [127] "GALNT6" "GPR65" "BCL2A1" "GOLGA6D" "ITGAX" "PIK3R5"
## [133] "VAV1" "MYO1F" "CD53" "FCGR2A" "LEFTY2" "PTPN7"
## [139] "OSBPL10" "TM4SF19" "PLA2G7" "DOK2" "CDKN2B" "GSN"
## [145] "RGS10" "ST14" "TAGLN" "FERMT3" "FCGR1A" "BMP6"
## [151] "SAMSN1" "SLC7A7" "VSIG4" "TDRD9" "CD109" "SLC34A2"
## [157] "TNFRSF14" "NCF1" "SLAMF8" "FCER1G" "C1QC" "RUNX1"
## [163] "GOLGA6A" "ITGB2" "IL6R" "LAPTM5" "FCMR" "HPGDS"
## [169] "CTSS" "S100A11" "S100A9" "CCKAR" "FBLN2" "CCR1"
## [175] "TMEM144" "TAGAP" "SYK" "ALDH1A1" "CYBB" "HTRA1"
## [181] "SPINT1" "PLD4" "MFAP4" "C15orf48" "MS4A7" "MS4A14"
## [187] "LAIR1" "LDLRAD4" "CXCL8" "CD52" "ITGAM" "FPR1"
## [193] "GPR34" "C3AR1" "SYNPO" "KLHL6" "C1QB" "C1QA"
## [199] "OLR1" "TLR10" "TLR6" "LRRC25" "NUPR1" "MSC"
## [205] "SSC5D" "HCLS1" "S1PR5" "CLDN7" "UBA7" "BGN"
## [211] "CSF1R" "TMEM119" "APOBR" "EVI2B" "CD300LF" "FFAR4"
## [217] "ARHGAP30" "TRPV2" "AMTN" "DPYD" "TLR7" "SPN"
## [223] "ADAMTSL2" "RCSD1" "SUCNR1" "FCGR3A" "LST1" "STK38L"
## [229] "APOC2" "HLA-DMB" "CLEC5A" "CLEC19A" "PECAM1" "RAB7B"
## [235] "CCL3" "ADORA3"print('Correlation of IFNy results between MSN38 & WTC11')## [1] "Correlation of IFNy results between MSN38 & WTC11"cor(line.res$DEG_COiMg[['MSN38']][rownames(line.res$DEG_COiMg[['MSN38']]) %in% line.res$DEG_COiMg$overlap_msn38Xwtc11,]$log2FoldChange,
line.res$DEG_COiMg[['WTC11']][rownames(line.res$DEG_COiMg[['WTC11']]) %in% line.res$DEG_COiMg$overlap_msn38Xwtc11,]$log2FoldChange)## [1] 0.8945231LinexCOiMG
#LinexCOiMG
for (i in c('FOXJ1',
'CLDN5',
'IDO1',
'GBP5',
'MMP10',
'SOX2',
'RBFOX3')){
p <- ggplot(pca.df_upgraded,aes_string(x='LinexCOiMG',y=i, fill='LinexCOiMG')) +
#stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title=paste0("LinexCOiMG ", i, " expression"),x="LinexCOiMG", y = paste0(i,' expression')) +
theme_classic()
plot(p)
}
Interaction analysis
Interaction analysis of co-culture*IFNy
What happens if you co-culture AND stimulate with IFNy?
Variance analysis
pca_cor3_mod1$interaction <- as.factor(pca_cor3_mod1$COiMg):as.factor(pca_cor3_mod1$condition)
PCAplot(pca_cor3_mod1, "interaction", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2,
colors = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"Spectral")))(6)) 
DEA
Volcano
results$COiMGxIFNg$symbol <- rownames(results$COiMGxIFNg)
volcano_plot(data.frame(results$COiMGxIFNg), title = 'CTR vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 10, ymax2 = 11)
DEA boxplots
COiMGxCondition
for (i in c('FOXJ1',
'CLDN5',
'IDO1',
'GBP5',
'MMP10',
'SOX2',
'RBFOX3')){
p <- ggplot(pca.df_upgraded,aes_string(x='COiMGxCondition',y=i, fill='COiMGxCondition')) +
#stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title=paste0("COiMGxCondition ", i, " expression"),x="COiMGxCondition", y = paste0(i,' expression')) +
theme_classic()
plot(p)
}
Full LinexConditionxCOiMG boxplots
for (i in c('FOXJ1',
'CLDN5',
'IDO1',
'GBP5',
'MMP10',
'SOX2',
'RBFOX3')){
p <- ggplot(pca.df_upgraded_sub,aes_string(x='LinexConditionxCOiMg',y=i, fill='LinexConditionxCOiMg')) +
#stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
labs(title=paste0("LinexConditionxCOiMg ", i, " expression"),x="LinexConditionxCOiMg", y = paste0(i,' expression')) +
theme_classic()+
rotate_x_text(angle=45) +
theme_classic()
plot(p)
}
